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Abstract — An algorithm that performs joint equalization and 
decoding for channels with nonlinear two-dimensional inter- 
symbol interference is presented. The algorithm performs sum- 
product message-passing on a factor graph that represents the 
underlying system. The two-dimensional optical storage (Two- 
DOS) technology is an example of a system with nonlinear 
two-dimensional intersymbol interference. Simulations for the 
nonlinear channel model of TwoDOS show significant improve- 
ment in performance over uncoded performance. Noise tolerance 
thresholds for the TwoDOS channel computed using density 
evolution are also presented. 

I. Introduction 

Two-dimensional (2D) intersymbol interference (ISI) chan- 
nels have received a lot of attention lately. This is mainly 
due to the fact that research focus in the storage industry is 
shifting towards developing a two-dimensional paradigm for 
storage. Current storage technologies are restricted by physical 
limits which will prevent them from keeping up with the 
ever increasing demands for data storage. This has prompted 
the development of technologies that employ novel tech- 
niques for data storage. Patterned magnetic media, in which 
information is stored in isolated single grain islands, make 
areal densities of the order of 10 12 bits/in 2 feasible, which 
is far beyond the limits of conventional magnetic recording. 
Holographic storage, page-oriented optical memories, and the 
two-dimensional optical storage (TwoDOS) technology are 
potential optical storage technologies of the future. Due to 
the two-dimensional nature of storage, these advanced storage 
technologies have 2D ISI during the readback process. Con- 
ventional recording technologies like magnetic hard disks and 
DVDs have one-dimensional ISI for which partial response 
maximum-likelihood decoding has been highly successful. Ex- 
tending PRML to two dimensions is not straightforward since 
maximum-likelihood decoding in two-dimensions is computa- 
tionally infeasible. This motivates the need for new methods to 
combat 2D ISI. Besides advanced storage technologies, multi- 
user communication scenarios, like cellular communication, 
also have situations where 2D ISI is prevalent. 

Detection schemes for 2D ISI channels have been proposed 
by many researchers [l]-[5]. Singla et ai, [6], [7] have 
proposed joint equalization and decoding schemes for 2D ISI 
channels and have shown the benefit of using error-correction 
coding in conjunction with detection. More often than not, the 
ISI is modeled as a linear filter. Although a good starting point, 



the linearity assumption doesn't hold in general. TwoDOS is 
an example of a system where the ISI is nonlinear. 

TwoDOS is, potentially, the next generation optical storage 
technology with projected storage capacity twice that of the 
blu-ray disk and with ten times faster data access rates [8], 
[9]. As in conventional optical disk recording, bits in the 
TwoDOS model are written on the disk in spiral tracks. 
However, instead of having a single row of bit cells, each 
track consists of a number of bit rows stacked together. 
Thus, TwoDOS is a truly two-dimensional storage paradigm. 
Successive tracks on the disk are separated by a guard band 
which consists of one empty bit row. In addition, the bit cells 
are hexagonal; this allows 15 percent higher packing density 
than rectangular bit cells leading to higher storage capacity. 
As in conventional optical disk recording, a 0/1 is represented 
by the absence/presence of a pit on the disk surface. A scalar 
diffraction model proposed by Coene [8] for optical recording 
is used to model the readback signal from the disk. Under 
this model the readback intensity from the disk has linear and 
bilinear contributions from the stored data. 

Various detection schemes for TwoDOS have been pro- 
posed [9]-[ 1 1] . These schemes, with the exception of that 
proposed by Immink et ai, [9], use two-dimensional partial 
response equalization to obtain a linear channel model for 
the ISI. Then, equalization methods like minimum mean- 
squared-error equalization, are used for detection on this linear 
channel model. Since partial response equalization leads to 
noise correlation there is an inherent loss associated with these 
schemes. Thus, it is prudent to search for decoding schemes 
that avoid partial response equalization and are designed taking 
into account the nonlinear structure of the ISI. Immink et 
ah, [9] propose using a stripe-wise Viterbi detector that is 
designed for the nonlinear ISI. However, they did not employ 
any error correction coding. 

In this paper, a low-complexity scheme for joint equalization 
and decoding for nonlinear 2D ISI channels is presented. The 
scheme was first proposed for linear 2D ISI channels [6] and 
has been appropriately modified for the nonlinear channel. 
This scheme, called the full graph scheme, performs sum- 
product message-passing on a joint graph that represents the 
error-correction code and the nonlinear 2D ISI channel. Low- 
density parity-check (LDPC) codes [12] are used for error 
correction. Simulations for the nonlinear channel model of 



TwoDOS demonstrate the potential of using the full graph 
scheme. Significant improvement in performance is observed 
over uncoded performance. Noise tolerance thresholds are 
calculated for regular LDPC codes of different rates and full 
graph decoding for the nonlinear 2D ISI channel. 

The paper is organized as follows. The model of the 
system and the channel model for TwoDOS is described in 
Section lU The full graph message-passing algorithm and its 
performance for TwoDOS are presented in Section [HI] The 
density evolution algorithm and the noise tolerance thresholds 
are presented in Section [Tvl Section IV! concludes the paper. 

II. System model 

The system is modeled as a discrete-time communication 
system governed by the following equation; 

r(i,j) = h({x(k,l) : {k,l)eKj})+w(i,j). (1) 

Here, r(i, j) is the data received at the output of the channel. 
x(k, I) are the channel inputs, obtained by encoding the user 
data using an error correction code. The user data and the 
encoded data are assumed to be binary. LDPC codes are used 
for error correction. w(i,j) are samples of additive white 
Gaussian noise (AWGN) with zero mean and variance a 2 . Afij 
is the set of indices of all the bits that interfere with x(i,j) 
during readback and h(-) is the function that encapsulates the 
nonlinear 2D interference. 

For TwoDOS, a scalar diffraction model proposed by Co- 
ene [8] for optical recording is used to model the readback 
signal. Using the model, the readback signal (optical intensity) 
from the disk is 

r(i,j) = l-^2c(k,l)x(k,l) + 
(k,i) 

+ d(k,l;m,n)x(k,l)x(m,n), (2) 

(fc,Z)#(m,n) 

where c(k, I) and d(k, I; m, n) are respectively the linear and 
nonlinear ISI coefficients. These coefficients depend on the 
parameters of the optical system, such as the wavelength 
of the laser, numerical aperture of the readback lens and 
geometry of the recording (pit and track dimensions). The 
extent of the interference is limited by the spot size of the 
laser used for reading. For low-to-moderate storage densities 
this leads to interference from the nearest neighbors only, 
whereas, for high storage densities (twice that of the blu- 
ray disk) the interference from bits in the other 'shells' also 
becomes significant [8]. 

Using a nearest neighbor interference model, the signal 
intensity in depends on the data bit stored in the central 
bit cell and the 6 neighboring bit cells. If it is assumed 
that two configurations with the same central bit and same 
number of nonzero neighbors have identical signal values 
then the signal intensity takes on 14 values corresponding to 
the 14 different configurations. Fig. shows four of these 
14 configurations. As shown by Coene [8], this symmetry 



assumption is a good approximation. Table H] lists the signal 
levels for the 14 different configurations. 




Fig. 1. Four of the possible 14 nearest neighbor configurations for TwoDOS. 
The dark circles in the cells depict a pit corresponding to a stored 1 . Absence 
of a pit indicates a stored 0. The pits cover only about half the area of the 
hexagonal bit cells. This is done to reduce signal folding [8]. 



TABLE I 

Signal levels for TwoDOS recording using scalar diffraction 
model and nearest neighbors. (reproduced from [8]) 



Nonzero 


Central 


Central 


neighbors (n) 


bit=0 (s„ ) 


bit=l (s ni ) 





0.95 


0.50 


1 


0.80 


0.35 


2 


0.70 


0.30 


3 


0.55 


0.20 


4 


0.45 


0.15 


5 


0.35 


0.10 


6 


0.25 


0.05 



This table is reproduced from [8] and corresponds to 
TwoDOS recording with hexagonal lattice parameter and pit- 
hole diameter equal to 165 nm and 120 nm, respectively. 
Using these parameters gives TwoDOS a 1.4 fold increase 
in storage density over the blu-ray disc. Looking at the table, 
the nonlinearity of the ISI is quite apparent; the signal level 
does not change linearly as the number of nonzero neighbors 
increases and the range when the central bit is a is greater 
than when the central bit is 1. Similar to Immink et al., [9] 
the signal levels in Table H] take into account the interference 
from the outer shells (beyond the nearest neighbors) by taking 
an average across all the bit-patterns in these shells. 

III. Full Graph Message-Passing Algorithm 

The full graph algorithm computes approximate APPs of the 
codeword bits given the observations by performing message- 
passing on a factor graph representing the underlying system. 
This "full graph" has three types of nodes; variable nodes, 
check nodes, and measured data nodes corresponding to the 
codeword bits, the parity-check equations, and the observed 
data symbols, respectively. The upper two levels in the full 
graph represent the LDPC code bipartite graph showing how 
the codeword bits are connected to the check nodes via 
the LDPC code parity-check matrix. The lower two levels 
represent the channel ISI graph showing how the ISI induces 
dependencies between the codeword bits. Fig. |2] shows an il- 
lustration of the full graph where nearest neighbor interference 
is assumed. 

Message-passing on this full graph is performed using the 
following schedule: variable nodes to check nodes, check 
nodes to variable nodes, variable nodes to measured data nodes 
and finally measured data nodes to variable nodes. Following is 
a brief description of how the messages are updated for each of 




Fig. 2. Factor graph representation of the TwoDOS system. 

the aforementioned four steps for the sum-product algorithm. 
The messages in the update equations are probabilities. 
Variable-to-check messages: The message from a variable 
node a; to a check node c at the Zth iteration is calculated 
using the messages passed to x from its neighboring check 
and measured data nodes at the (I — l)th iteration. 

r£N r (x) c'£N c (x)\c 

/4ic(i) = « n n Mc^2(i),(3) 

re.N r (x) c'GiV c (3:)\c 

where /iH c (-)> Mrla:(')> an d £tclx(") are respectively, the 
variable-to-check, measured data-to-variable, and check-to- 
variable messages at the Zth iteration. N c (x) and N r (x) are the 
neighboring check and measured data nodes of x respectively, 
and a is a normalizing constant. 

Check-to-variable messages: At the Zth iteration the check- 
to-variable messages are calculated using the variable-to-check 
messages at the Zth iteration and the sum-product rule, 

/4i x (0) = J2 P(o\x = 0,{x'})- 

{x':x'£N(c)\x} 

• n 

x'eN(c)\x 

MEx(l) = E P(c\x = l,W}> 

{x':x'&N(c)\x} 

• n ^Lc(^), (4) 

x'£N(c)\x 

where N(c) are the variable nodes connected to check node 
c. 

Variable-to-measured data messages: The variable-to- 
measured data messages at the Zth iteration are calculated using 
the messages received at the variable nodes from the check 
nodes at the Zth iteration and from the measured data nodes 
at the (Z — l)th iteration, 



/4i r (o) = a n M^uo) n 4^(0) 

c£N c (x) r'eN r (x)\r 

,(0 m _ „ TT „(0 TT „(*-!)/ 



41,(1) 



« n ^-ix(i) n (5) 

c£JV c (i) r'eJV r (x)\r 

Measured data-to-variable messages: To complete one itera- 
tion, the measured data-to-variable messages at the Zth iteration 
are computed using the variable-to-measured data messages at 
the Zth iteration and the sum-product rule, 

,41.(0) = a £ p(r\x = 0,{x'}y 

{x':x'eN(r)\x} 

■ II ^ri*') 
x'£N(r)\x 

^l x (l) = a E p(r| a! = l J {x'})- 

{a/:x'eiV(r)\x} 

• n MxU(^). (6) 

i'eA r (r)\i 

N(r) is the set of all variable nodes connected to measured 
data node r. The conditional probabilities above are values 
of a Gaussian probability density function. The mean of this 
probability density function is determined by the signal levels 
given in Table|I]and the variance is equal to the noise variance. 

After this step the "pseudo-posterior" probabilities are cal- 
culated using the messages from the check nodes and the 
measured data nodes, 

«d°(o) = a n m&.(o) n ^cix(o) 

r£N r (x) c£N c (x) 

^(i) = a n ^x(i) n ^x(i)- (7) 

r£N r (x) c£N a (x) 

The codeword estimate is obtained by setting x to 1 if 
<j4^(l) > qx\o) and otherwise. The decoding stops if the 
decoder converges to a codeword or a maximum number of 
iterations is exhausted. The complexity of the full graph algo- 
rithm is linear in the LDPC code block length and quadratic 
in the size of the interfering neighborhood. 

Results of using full graph message-passing for the Two- 
DOS channel model of (0 and Table [I] are shown in Fig. [3] 
The LDPC code used is a block length 10000, regular (3,30) 
code. This high-rate code is chosen so as to add only a small 
number of redundant parity bits. The SNR (E^/No) is defined 
as the average signal energy divided by the noise power; 



SNR = lOlogi 



En=0 (n) ( S 



S 2 
n 



/2 7 



R-(2a 2 ) 



(8) 



where s n Js ni is the signal level given that the central bit is 
a 0/1 and has n nonzero neighbors and R is the LDPC code 
rate. From right to left the curves in Fig. [3] are; the uncoded 
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Fig. 3. Results of using the full graph message-passing algorithm for the 
TwoDOS recording model given by |2) and Table Q 



performance, the performance of the full graph message- 
passing algorithm for a maximum of 1, 2, 3, 4, and 5 iterations. 
The number of iterations is kept low so as to reduce decoding 
delay. As the curves show, the improvement in performance is 
quite significant. The performance of the full graph algorithm 
after 5 iterations is about 8 dB better (at a bit error rate of 
10 -6 ) than the performance reported in [9], where a stripe- 
wise Viterbi detector is used for the same ISI. 

For the results shown here the full graph algorithm only 
takes into account the ISI from the nearest neighbors. Includ- 
ing the interference from the other shells, though computa- 
tionally expensive, can still be accomplished at a significantly 
lower cost than the stripe-wise Viterbi detector where the 
complexity increases exponentially. Immink et al., [9] have 
proposed the use of modulation coding to eliminate data 
patterns that severely degrade the performance of the stripe- 
wise Viterbi detector. Furthermore, they have incorporated 
the modulation code and detection into a Bliss-like scheme. 
The modulation code can be incorporated seamlessly into 
our system model; the modulation code can be represented 
graphically and will add another level to the full graph. The 
connections of this subgraph will differ depending on whether 
the modulation code is incorporated in the usual way (after 
the LDPC encoder) or used in a Bliss-like scheme. We have 
not yet studied the performance of such a system. 

IV. Density Evolution and Threshold 
Computation 

For many channels and decoders of interest, LDPC codes 
exhibit a threshold phenomenon [12]; there exists a critical 
value of the channel parameter (noise tolerance threshold, 
say 6*) such that an arbitrarily small bit error probability 
can be achieved if the noise level (r5) is smaller than 6* 
and the code length is long enough. On the other hand, for 
5 > 5* the probability of bit error is larger than a positive 



constant. Richardson and Urbanke [12] developed an algorithm 
called density evolution for iteratively calculating message 
densities, enabling the determination of the aforementioned 
threshold. Kavcic et al, [13] extended the work of Richardson 
and Urbanke to compute noise tolerance thresholds for one- 
dimensional ISI channels. 

Using a density evolution similar to that proposed by Kavcic 
et al, [13] noise tolerance thresholds are computed for the 
full graph algorithm and the nonlinear TwoDOS channel. 
Following is a brief description of the algorithm. Density 
evolution tracks the "evolution" of the probability density 
function (pdf) of correct (or incorrect) messages passed on the 
graph as the iterations progress. Density evolution is described 
more conveniently using log-likelihood ratios (LLR) instead 
of probabilities to represent the messages. The LLR for the 
variable-to-check messages at the Zth iteration is defined as 



-tvx— >c — lOg 77 



= (0) 



The LLRs L r 



and L r 



are defined analogously. Using LLRs gives an equivalent 
representation of the full graph message-passing algorithm. 
Equations Q-© can be rewritten as 



L (l) 



E L l 

r£N r (x) 



E i 

c'£N c (x)\c 



2(9) 



tanh 



L 



(0 



2 

L (i)_ 



= (-l) c J] tanh 



r(0 

x'^c I (1Q) 



x'EN{c)\x 

= E 41* + E J&fri) 

c£N c (x) r'eN r (x)\r 

41 x = F(L^ r ;x'eN(r)\x), (12) 

The update for measured data-to-variable messages has no 
closed form when represented using LLRs. The "tanh" 
rule [14] cannot be applied to since the measured data 
nodes are not binary-valued. Equation H2\ represents the 
measured data-to-variable message update via a function F 
which performs the appropriate computation. 

Let f$p (x) be the pdf of the correct message from a variable 
node to a check node at the Zth round of message-passing. This 
pdf is evolved through ll9l-J12>. The evolution of the density 
functions through l|9jl and i ll 11 are simple convolutions and can 
be implemented efficiently using the fast Fourier transform. 
Density evolution for J 1 Oft can be implemented using the 
change in measure described in [12] or more efficiently by 
using a table-lookup as explained in [15]. For density evolution 
through fl!2i Monte Carlo simulations are used; message- 
passing is performed on the channel graph using a long 
block length and the pdf of the outgoing messages from the 
measured data nodes is approximated by using the histogram 
of computed messages. 

After evolving fP(x) through §-0, fi l+1) (x) is ob- 
tained. Using fv (x), the error probability at the (I + l)th 
iteration, pi 1 +1 \ can be computed as 



fi l+1) {x)dx. 



(13) 



The noise tolerance threshold, 5*, can be calculated as the 
supremum of all S for which the error probability goes to 
zero as the iterations progress. Table |n] shows the computed 
thresholds for regular LDPC codes of different rates. The noise 
tolerance threshold is the variance of the AWGN. The SNR 
is calculated as in l|8} except that the rate of the code is not 
taken into account. It should be noted that to compensate for 
the loss in density due to error-correction coding the bits need 
to be stacked closer, leading to increased ISI. This effect is 
not taken into account during the calculation of the thresholds 
for the different rate codes. 

TABLE II 

Noise tolerance thresholds for the TwoDOS channel. 



LDPC 


Code 


Threshold 


Threshold 


Code 


Rate 


°i 


SNR [dB] 


(3,3) 


0.000 


0.0670 


1.6070 


(3,4) 


0.250 


0.0436 


3.4729 


(3,5) 


0.400 


0.0283 


5.3499 


(3,6) 


0.500 


0.0215 


6.5433 


(3,9) 


0.667 


0.0140 


8.4064 


(3,12) 


0.750 


0.0117 


9.1859 


(3,15) 


0.800 


0.0103 


9.7393 


(3,30) 


0.900 


0.0061 


12.0144 


(3,60) 


0.950 


0.0035 


14.4270 


(3,90) 


0.967 


0.0030 


15.0965 


(3,120) 


0.975 


0.0027 


15.5541 


(3,150) 


0.980 


0.0025 


15.8883 



In proving the existence of thresholds for memoryless chan- 
nels the crucial innovation of Richardson and Urbanke [12] 
was the "concentration results." These results state that as the 
block length tends to infinity the performance of the LDPC 
decoder on random graphs converges to its expected behavior 
and that the expected behavior can be determined from the 
corresponding cycle-free behavior. Kavcic et al., [13] extended 
these concentration results to one-dimensional ISI channels by 
using LDPC coset codes. 

For 2D ISI channels the concentration results do not hold 
since the channel graph has short cycles even in the limit 
of infinitely long block length. Hence existence of thresholds 
cannot be proved using the concentration analysis. However, 
our simulations seem to suggest that for the TwoDOS channel 
the full graph algorithm respects the thresholds computed 
using density evolution. Simulations for long block lengths 
(^2xl0 6 ) show that very low bit error rates (^10~ 6 ) are 
obtained only when the noise variance is smaller than the 
threshold. Although this does not prove the existence of a 
threshold, it suggests that the noise tolerance thresholds of 
Table ITT1 are upper bounds on the performance of the full graph 
algorithm. Besides that, the thresholds also serve as a design 
parameter; given a system with a specified SNR it is sufficient 
to pick an LDPC code having a smaller threshold SNR thereby 
ensuring that the bit-error rate can be made arbitrarily small 
as the block length increases. 

V. Conclusion 

A message-passing based scheme for joint equalization and 
decoding for nonlinear two-dimensional intersymbol interfer- 



ence channels has been proposed. The scheme, called the full 
graph algorithm, performs sum-product message-passing on a 
joint graph of the error correction code and the channel. The 
complexity of the full graph algorithm is linear in the block 
length of the error correction code and quadratic in the size of 
interference neighborhood. The performance of the algorithm 
is studied for the two-dimensional optical storage paradigm. 
Simulations for the nonlinear channel model of TwoDOS 
show significant improvement over uncoded performance. The 
performance is about 8 dB better than that reported by Immink 
et al, [9] for the same ISI. Using density evolution noise 
tolerance thresholds for the full graph algorithm are also 
computed. 
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